      program main

!*********************************************************************72
!
!! MAIN is the main program for RANLIB_PRB.
!
!  Licensing:
!
!    This code is distributed under the GNU LGPL license.
!
!  Modified:
!
!    02 April 2013
!
!  Author:
!
!    Original FORTRAN77 version by Barry Brown, James Lovato.
!    Modifications by John Burkardt
!
      implicit none

	real a, r, r2, mu, sigmasq, irand, aa(4), rr(4) 
	real MeanIn(4), VarIn(4), gengam, array(2000000)
	real CumProb(4), av, var, xmin, xmax, sgamma
      integer*4 i, j, m, n 
      character * ( 100 ) phrase

      call timestamp ( )
      write ( *, '(a)' ) ' '
      write ( *, '(a)' ) 'RANLIB_PRB'
      write ( *, '(a)' ) '  FORTRAN77 version'
      write ( *, '(a)' ) '  Test the RANLIB library.'

      CALL init_random_seed()         ! see example of RANDOM_SEED

      CALL RANDOM_NUMBER(r2)

OPEN(UNIT=12, FILE="GammaRNsout.txt", ACTION="write", STATUS="replace")
OPEN(UNIT=13, FILE="GammaRNsout2.txt", ACTION="write", STATUS="replace")

!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
!
!  Four levels of claim frequency/severity are parametrized
!   I.   High frequency low severity MeanIn(1); VarIn(1)
!   II.  Moderate frequency moderate severity MeanIn(2); VarIn(2)
!   I.   Low frequency high severity MeanIn(3); VarIn(3)
!   I.   Rare, catastrophic severity MeanIn(4); VarIn(4)


MeanIn(1) = 50.
MeanIn(2) = 200.
MeanIn(3) = 25000.
MeanIn(4) = 500000.
VarIn(1)  = 10000000.
VarIn(2)  = 10000000.
VarIn(3)  = 50000000.
VarIn(4)  = 80000000.


!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
!
!	Gamma Parameters for four types of claims calculated from 
!	mean and variance of the claim frequency/severity distribution
!
!
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!

do i = 1,4
	aa(i) = MeanIn(i)/VarIn(i)
	rr(i) = MeanIn(i)**2/VarIn(i)
	write (*, '(a, 4g12.6)') 'MeanIn(1), VarIn(1), aa, rr  ', MeanIn(i), VarIn(i), aa(i), rr(i)
enddo !
!
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!

!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
!
!  Four levels of claim frequency probability
!   I.   High frequency low severity :: Cumulative Probability 0.45000
!   II.  Moderate frequency moderate severity :: Cumulative Probability 0.99000
!   I.   Low frequency high severity :: Cumulative Probability 0.99500
!   I.   Rare, catastrophic severity :: Cumulative Probability 1.00000
!
!   

   CumProb(1) = 0.35000
   CumProb(2) = 0.98000
   CumProb(3) = 0.99500
   CumProb(4) = 1.00000
!
!
!
!
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!

!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
!
!
!	Generate a sample of 'n' claims from the four tiers of claims
!
!
m = 1000

do j = 1, m
	write (*, '(a, i8, 4f12.0)') 'First Loop Number m: ', j
	n = 100000
	do i = 1, n
		write (*, '(a, i8, 4f12.0)') 'Second Loop Number n: ', i
	     CALL RANDOM_NUMBER(irand)

		IF (irand .lt. CumProb(1)) THEN
			a = aa(1)
			r = rr(1)
		ELSE IF (irand .lt. CumProb(2)) THEN
			a = aa(2)
			r = rr(2)
		ELSE IF (irand .lt. CumProb(3)) THEN
			a = aa(3)
			r = rr(3)
		ELSE 
			a = aa(4)
			r = rr(4)
		END IF
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
	write (*, '(a, 2i8, 4f12.4)') 'First/Second Loop Numbers j, i, irand: ', j, i, a, r, irand

	array(i) = sgamma ( r ) / a

	write (12, '(f12.0)') array(i)

	enddo

      	call stats ( array, n, av, var, xmin, xmax )
	write (*, '(a, i8, 4f12.0)') 'Claim Number n, av, var, xmin, xmax: ', n, av, var, xmin, xmax

enddo

!
!  Terminate.
!
      write ( *, '(a)' ) ' '
      write ( *, '(a)' ) 'RANLIB_PRB'
      write ( *, '(a)' ) '  Normal end of execution.'
      write ( *, '(a)' ) ' '
      call timestamp ( )
close (12)
close (13)
      stop
      end


      SUBROUTINE init_random_seed()
        INTEGER :: i, n, clock
        INTEGER, DIMENSION(:), ALLOCATABLE :: seed
        
        CALL RANDOM_SEED(size = n)
        ALLOCATE(seed(n))
          
        CALL SYSTEM_CLOCK(COUNT=clock)
          
          seed = clock + 37 * (/ (i - 1, i = 1, n) /)
        CALL RANDOM_SEED(PUT = seed)
          
        DEALLOCATE(seed)
      END SUBROUTINE




subroutine timestamp ( )

!*****************************************************************************80
!
!! TIMESTAMP prints the current YMDHMS date as a time stamp.
!
!  Example:
!
!    May 31 2001   9:45:54.872 AM
!
!  Licensing:
!
!    This code is distributed under the GNU LGPL license. 
!
!  Modified:
!
!    31 May 2001
!
!  Author:
!
!    John Burkardt
!
!  Parameters:
!
!    None
!
  implicit none

  character ( len = 8 ) ampm
  integer ( kind = 4 ) d
  character ( len = 8 ) date
  integer ( kind = 4 ) h
  integer ( kind = 4 ) m
  integer ( kind = 4 ) mm
  character ( len = 9 ), parameter, dimension(12) :: month = (/ &
    'January  ', 'February ', 'March    ', 'April    ', &
    'May      ', 'June     ', 'July     ', 'August   ', &
    'September', 'October  ', 'November ', 'December ' /)
  integer ( kind = 4 ) n
  integer ( kind = 4 ) s
  character ( len = 10 ) time
  integer ( kind = 4 ) values(8)
  integer ( kind = 4 ) y
  character ( len = 5 ) zone

  call date_and_time ( date, time, zone, values )

  y = values(1)
  m = values(2)
  d = values(3)
  h = values(5)
  n = values(6)
  s = values(7)
  mm = values(8)

  if ( h < 12 ) then
    ampm = 'AM'
  else if ( h == 12 ) then
    if ( n == 0 .and. s == 0 ) then
      ampm = 'Noon'
    else
      ampm = 'PM'
    end if
  else
    h = h - 12
    if ( h < 12 ) then
      ampm = 'PM'
    else if ( h == 12 ) then
      if ( n == 0 .and. s == 0 ) then
        ampm = 'Midnight'
      else
        ampm = 'AM'
      end if
    end if
  end if

  write ( *, '(a,1x,i2,1x,i4,2x,i2,a1,i2.2,a1,i2.2,a1,i3.3,1x,a)' ) &
    trim ( month(m) ), d, y, h, ':', n, ':', s, '.', mm, trim ( ampm )

  return
end


     function sgamma ( a )

!*********************************************************************72
!
!c SGAMMA samples the standard Gamma distribution.
!
!  Discussion:
!
!    This procedure corresponds to algorithm GD in the reference.     
!                                                                      
!  Licensing:
!
!    This code is distributed under the GNU LGPL license.
!
!  Modified:
!
!    01 April 2013
!
!  Author:
!
!    Original FORTRAN77 version by Barry Brown, James Lovato.
!!   Modifications by John Burkardt.
!
!  Reference:
!                                                                                                                                           
!    Joachim Ahrens, Ulrich Dieter,
!    Generating Gamma Variates by a Modified Rejection Technique,
!    Communications of the ACM,
!    Volume 25, Number 1, January 1982, pages 47-54.               
!1                                                                      
!1  Parameters:
!
!    Input, real A, the parameter of the standard gamma distribution.
!    0.0 < A < 1.0.
!
!    Output, real SGAMMA, a random deviate from the distribution.
!
      implicit none

      real a
      real a1
      real a2
      real a3
      real a4
      real a5
      real a6
      real a7
      real b
      real c
      real d
      real e
      real e1
      real e2
      real e3
      real e4
      real e5
      real p
      real q
      real q0
      real q1
      real q2
      real q3
      real q4
      real q5
      real q6
      real q7
      real r
      real r2
      real s
      real s2
      real sexpo
      real si
      real sgamma
      real snorm
      real sqrt32
      real t
      real u
      real v
      real w
      real x

      save a1
      save a2
      save a3
      save a4
      save a5
      save a6
      save a7
      save e1
      save e2
      save e3
      save e4
      save e5
      save q1
      save q2
      save q3
      save q4
      save q5
      save q6
      save q7
      save sqrt32

      data a1 /  0.3333333E+00 /
      data a2 / -0.2500030E+00 /
      data a3 /  0.2000062E+00 /
      data a4 / -0.1662921E+00 /
      data a5 /  0.1423657E+00 /
      data a6 / -0.1367177E+00 /
      data a7 /  0.1233795E+00 /
      data e1 / 1.0E+00 /
      data e2 / 0.4999897E+00 /
      data e3 / 0.1668290E+00 /
      data e4 / 0.0407753E+00 /
      data e5 / 0.0102930E+00 /
      data q1 /  0.04166669E+00 /
      data q2 /  0.02083148E+00 /
      data q3 /  0.00801191E+00 /
      data q4 /  0.00144121E+00 /
      data q5 / -0.00007388E+00 /
      data q6 /  0.00024511E+00 /
      data q7 /  0.00024240E+00 /
      data sqrt32 / 5.656854E+00 /
!
!  Recalculations if A has changed.
!
      if ( 1.0E+00 .le. a ) then

        s2 = a - 0.5E+00
        s = sqrt ( s2 )
        d = sqrt32 - 12.0E+00 * s
!
!  Immediate acceptance.
!
        t = snorm ( )
        x = s + 0.5E+00 * t
        sgamma = x * x

        if ( 0.0E+00 .le. t ) then
          return
        end if
!!
!  Squeeze acceptance.
!
        CALL RANDOM_NUMBER(u)
        if ( d * u .le. t * t * t ) then
          return
        end if
!
!  Possible recalculations.
!
        r = 1.0E+00 / a
        q0 = (((((( q7 * r + q6 ) * r + q5 ) * r + q4 ) * r + q3 ) * r + q2 ) * r + q1 ) * r
!
!  Approximation depending on size of parameter A.
!
        if ( 13.022E+00 .lt. a ) then
          b = 1.77E+00
          si = 0.75E+00
          c = 0.1515E+00 / s
        else if ( 3.686E+00 .lt. a ) then
          b = 1.654E+00 + 0.0076E+00 * s2
          si = 1.68E+00 / s + 0.275E+00
          c = 0.062E+00 / s + 0.024E+00
        else
          b = 0.463E+00 + s + 0.178E+00 * s2
          si = 1.235E+00
          c = 0.195E+00 / s - 0.079E+00 + 0.16E+00 * s
        end if
!
!  Quotient test.
!
        if ( 0.0E+00 .lt. x ) then

          v = 0.5E+00 * t / s

          if ( 0.25E+00 .lt. abs ( v ) ) then
            q = q0 - s * t + 0.25E+00 * t * t + 2.0E+00 * s2 * alog ( 1.0E+00 + v )
          else
            q = q0 + 0.5E+00 * t * t * (((((( a7 * v + a6 ) * v + a5 ) * v + a4 ) * v + a3 ) * v + a2 ) * v + a1 ) * v
          end if

          if ( alog ( 1.0E+00 - u ) .le. q ) then
            return
          end if

        end if

10      continue

        e = sexpo ( )
CALL RANDOM_NUMBER(r)
        u = 2.0E+00 * r - 1.0E+00

        if ( 0.0E+00 .le. u ) then
          t = b + abs ( si * e )
        else
          t = b - abs ( si * e )
        end if
!
!  Possible rejection.
!
        if ( t .lt. -0.7187449E+00 ) then
          go to 10
        end if
!
!  Calculate V and quotient Q.
!
        v = 0.5E+00 * t / s

        if ( 0.25E+00 .lt. abs ( v ) ) then
          q = q0 - s * t + 0.25E+00 * t * t + 2.0E+00 * s2 * alog ( 1.0E+00 + v )
        else
          q = q0 + 0.5E+00 * t * t * (((((( a7 * v + a6 ) * v + a5 ) * v + a4 )* v + a3 )* v + a2 ) * v + a1 ) * v
        end if
!
!  Hat acceptance.
!
        if ( q .le. 0.0E+00 ) then
          go to 10
        end if

        if ( 0.5E+00 .lt. q ) then
          w = exp ( q ) - 1.0E+00
        else
          w = (((( e5 * q + e4 ) * q + e3 ) * q + e2 ) * q + e1 ) * q
        end if
!
!  May have to sample again.
!
        if ( w * exp ( e - 0.5E+00 * t * t ) .lt. c * abs ( u ) ) then
          go to 10
        end if

        x = s + 0.5E+00 * t
        sgamma = x * x
        return
!
!  Method for A < 1.
!
      else

        b = 1.0E+00 + 0.3678794E+00 * a

20      continue

CALL RANDOM_NUMBER(r)
        p = b * r

        if ( p .lt. 1.0E+00 ) then

          sgamma = exp ( alog ( p ) / a )

          if ( sgamma .le. sexpo ( ) ) then
            return
          end if

          go to 20

        end if

        sgamma = - alog ( ( b - p ) / a )

        if ( sexpo ( ) .lt. ( 1.0E+00 - a ) * alog ( sgamma ) ) then
          go to 20
        end if

      end if

      return
      end

      function genunf ( low, high )

!c*********************************************************************72
!c
!cc GENUNF generates a uniform random deviate.
!c
!c  Discussion:
!c
!c    This procedure generates a real deviate uniformly distributed between 
!1c    LOW and HIGH.
!c
!c  Licensing:
!c
!c    This code is distributed under the GNU LGPL license.
!c
!c  Modified:
!c
!c    20 March 2013
!c
!c  Author:
!c
!c    Original FORTRAN77 version by Barry Brown, James Lovato.
!c    Modifications by John Burkardt.
!c
!c  Parameters:
!c
!c    Input, real LOW, HIGH, the lower and upper bounds.
!c
!c    Output, real GENUNF, a random deviate from the distribution.
!c 
      implicit none

      real genunf
      real high
      real low
      real r

CALL RANDOM_NUMBER(r)
      genunf = low + ( high - low ) * r

      return
      end


subroutine stats ( x, n, av, var, xmin, xmax )

!*****************************************************************************80
!
!! STATS computes statistics for a given array.
!
!  Discussion:
!
!    This procedure computes the average and variance of an array.
!
!  Licensing:
!
!    This code is distributed under the GNU LGPL license.
!
!  Modified:
!
!    31 March 2013
!
!  Author:
!
!    Original FORTRAN77 version by Barry Brown, James Lovato.
!    FORTRAN90 version by John Burkardt.
!
!  Parameters:
!
!    Input, real ( kind = 4 ) X(N), the array to be analyzed.
!
!    Input, integer ( kind = 4 ) N, the dimension of the array.
!
!    Output, real ( kind = 4 ) AV, the average value.
!
!    Output, real ( kind = 4 ) VAR, the variance.
!
!    Output, real ( kind = 4 ) XMIN, XMAX, the minimum and maximum entries.
!
  implicit none

  integer ( kind = 4 ) n

  real ( kind = 4 ) av
  integer ( kind = 4 ) i
  real ( kind = 4 ) total
  real ( kind = 4 ) var
  real ( kind = 4 ) x(n)
  real ( kind = 4 ) xmax
  real ( kind = 4 ) xmin

  xmin = x(1)
  xmax = x(1)
  total = 0.0E+00
  do i = 1, n
    total = total + x(i)
    xmin = min ( xmin, x(i) )
    xmax = max ( xmax, x(i) )
  end do

  av = total / real ( n )
! write (*, '(a, i8, f12.0)') 'Number claims and total claims: ', n, total
  total = 0.0E+00
  do i = 1, n
    total = total + ( x(i) - av ) ** 2
  end do
  var = total / real ( n - 1 )

  return
end
subroutine trstat ( pdf, parin, av, var )

!*****************************************************************************80
!
!! TRSTAT returns the mean and variance for distributions.
!
!  Discussion:
!
!    This procedure returns the mean and variance for a number of statistical
!    distributions as a function of their parameters.
!
!    The input vector PARIN is used to pass in the parameters necessary
!    to specify the distribution.  The number of these parameters varies
!    per distribution, and it is necessary to specify an ordering for the
!    parameters used to a given distribution.  The ordering chosen here
!    is as follows:
!
!    bet
!      PARIN(1) is A
!      PARIN(2) is B
!    bin
!      PARIN(1) is Number of trials
!      PARIN(2) is Prob Event at Each Trial
!    chi
!      PARIN(1) = df
!    exp
!      PARIN(1) = mu
!    f
!      PARIN(1) is df numerator
!      PARIN(2) is df denominator
!    gam
!      PARIN(1) is A
!      PARIN(2) is R
!    nbn
!      PARIN(1) is N
!      PARIN(2) is P
!    nch
!      PARIN(1) is df
!      PARIN(2) is noncentrality parameter
!    nf
!      PARIN(1) is df numerator
!      PARIN(2) is df denominator
!      PARIN(3) is noncentrality parameter
!    nor
!      PARIN(1) is mean
!      PARIN(2) is standard deviation
!    poi
!      PARIN(1) is Mean
!    unf
!      PARIN(1) is LOW bound
!      PARIN(2) is HIGH bound
!
!  Licensing:
!
!    This code is distributed under the GNU LGPL license.
!
!  Modified:
!
!    01 April 2013
!
!  Author:
!
!    Original FORTRAN77 version by Barry Brown, James Lovato.
!    FORTRAN90 version by John Burkardt.
!
!  Parameters:
!
!    Input, character * ( 4 ) PDF, indicates the distribution:
!    'bet'  beta distribution
!    'bin'  binomial
!    'chi'  chisquare
!    'exp'  exponential
!    'f'    F (variance ratio)
!    'gam'  gamma
!    'nbn'  negative binomial
!    'nch'  noncentral chisquare
!    'nf'   noncentral f
!    'nor'  normal
!    'poi'  Poisson
!    'unf'  uniform
!
!    Input, real ( kind = 4 ) PARIN(*), the parameters of the distribution.
!
!    Output, real ( kind = 4 ) AV, the mean of the specified distribution.
!
!    Output, real ( kind = 4 ) VAR, the variance of the specified distribuion.
!
  implicit none

  real ( kind = 4 ) a
  real ( kind = 4 ) av
  real ( kind = 4 ) b
  integer ( kind = 4 ) n
  real ( kind = 4 ) p
  real ( kind = 4 ) parin(*)
  character * ( 4 ) pdf
  real ( kind = 4 ) r
  real ( kind = 4 ) var
  real ( kind = 4 ) width

  if ( pdf == 'bet' ) then

    av = parin(1) / ( parin(1) + parin(2) )
    var = ( av * parin(2) ) / ( ( parin(1) + parin(2) ) * &
      ( parin(1) + parin(2) + 1.0E+00 ) )

  else if ( pdf == 'bin' ) then

    n = int ( parin(1) )
    p = parin(2)
    av = real ( n ) * p
    var = real ( n ) * p * ( 1.0E+00 - p )

  else if ( pdf == 'chi' ) then

    av = parin(1)
    var = 2.0E+00 * parin(1)

  else if ( pdf == 'exp' ) then

    av = parin(1)
    var = av ** 2

  else if ( pdf == 'f' ) then

    if ( parin(2) <= 2.0001E+00 ) then
      av = -1.0E+00
    else
      av = parin(2) / ( parin(2) - 2.0E+00 )
    end if

    if ( parin(2) <= 4.0001E+00 ) then
      var = -1.0E+00
    else
      var = ( 2.0E+00 * parin(2) ** 2 * ( parin(1) + parin(2) - 2.0E+00 ) ) / &
        ( parin(1) * ( parin(2) - 2.0E+00 ) ** 2 * ( parin(2) - 4.0E+00 ) )
    end if

  else if ( pdf == 'gam' ) then

    a = parin(1)
    r = parin(2)
    av = r / a
    var = r / a ** 2

  else if ( pdf == 'nbn' ) then

    n = int ( parin(1) )
    p = parin(2)
    av = n * ( 1.0E+00 - p ) / p
    var = n * ( 1.0E+00 - p ) / p ** 2

  else if ( pdf == 'nch' ) then

    a = parin(1) + parin(2)
    b = parin(2) / a
    av = a
    var = 2.0E+00 * a * ( 1.0E+00 + b )

  else if ( pdf == 'nf' ) then

    if ( parin(2) <= 2.0001E+00 ) then
      av = -1.0E+00
    else
      av = ( parin(2) * ( parin(1) + parin(3) ) ) &
        / ( ( parin(2) - 2.0E+00 ) * parin(1) )
    end if

    if ( parin(2) <= 4.0001E+00 ) then
      var = -1.0E+00
    else
      a = ( parin(1) + parin(3) ) ** 2 &
        + ( parin(1) + 2.0E+00 * parin(3) ) * ( parin(2) - 2.0E+00 )
      b = ( parin(2) - 2.0E+00 ) ** 2 * ( parin(2) - 4.0E+00 )
      var = 2.0E+00 * ( parin(2) / parin(1) ) ** 2 * ( a / b )
    end if

  else if ( pdf == 'nor' ) then

    av = parin(1)
    var = parin(2) ** 2

  else if ( pdf == 'poi' ) then

    av = parin(1)
    var = parin(1)

  else if ( pdf == 'unf' ) then

    width = parin(2) - parin(1)
    av = parin(1) + width / 2.0E+00
    var = width ** 2 / 12.0E+00

  else

    write ( *, '(a)' ) ' '
    write ( *, '(a)' ) 'TRSTAT - Fatal error!'
    write ( *, '(a)' ) '  Illegal input value for PDF.'
    stop

  end if

  return
end

function snorm ( )

!*****************************************************************************80
!
!! SNORM samples the standard normal distribution.
!
!  Discussion:
!
!    This procedure corresponds to algorithm FL, with M = 5, in the reference.
!
!  Licensing:
!
!    This code is distributed under the GNU LGPL license.
!
!  Modified:
!
!    31 March 2013
!
!  Author:
!
!    Original FORTRAN77 version by Barry Brown, James Lovato.
!    FORTRAN90 version by John Burkardt.
!
!  Reference:
!
!    Joachim Ahrens, Ulrich Dieter,
!    Extensions of Forsythe's Method for Random
!    Sampling from the Normal Distribution,
!    Mathematics of Computation,
!    Volume 27, Number 124, October 1973, page 927-937.
!
!  Parameters:
!
!    Output, real ( kind = 4 ) SNORM, a random deviate from the distribution.
!
  implicit none

  real ( kind = 4 ) a(32)
  real ( kind = 4 ) aa
  real ( kind = 4 ) d(31)
  real ( kind = 4 ) h(31)
  integer ( kind = 4 ) i
  real ( kind = 4 ) r
  real ( kind = 4 ) s
  real ( kind = 4 ) snorm
  real ( kind = 4 ) t(31)
  real ( kind = 4 ) tt
  real ( kind = 4 ) u
  real ( kind = 4 ) ustar
  real ( kind = 4 ) w
  real ( kind = 4 ) y

  save a
  save d
  save h
  save t

  data a / &
        0.0000000E+00, 0.3917609E-01, 0.7841241E-01, 0.1177699E+00, &
        0.1573107E+00, 0.1970991E+00, 0.2372021E+00, 0.2776904E+00, &
        0.3186394E+00, 0.3601299E+00, 0.4022501E+00, 0.4450965E+00, &
        0.4887764E+00, 0.5334097E+00, 0.5791322E+00, 0.6260990E+00, &
        0.6744898E+00, 0.7245144E+00, 0.7764218E+00, 0.8305109E+00, &
        0.8871466E+00, 0.9467818E+00, 1.009990E+00,  1.077516E+00, &
        1.150349E+00,  1.229859E+00,  1.318011E+00,  1.417797E+00, &
        1.534121E+00,  1.675940E+00,  1.862732E+00,  2.153875E+00 /

  data d / &
        0.0000000E+00, 0.0000000E+00, 0.0000000E+00, 0.0000000E+00, &
        0.0000000E+00, 0.2636843E+00, 0.2425085E+00, 0.2255674E+00, &
        0.2116342E+00, 0.1999243E+00, 0.1899108E+00, 0.1812252E+00, &
        0.1736014E+00, 0.1668419E+00, 0.1607967E+00, 0.1553497E+00, &
        0.1504094E+00, 0.1459026E+00, 0.1417700E+00, 0.1379632E+00, &
        0.1344418E+00, 0.1311722E+00, 0.1281260E+00, 0.1252791E+00, &
        0.1226109E+00, 0.1201036E+00, 0.1177417E+00, 0.1155119E+00, &
        0.1134023E+00, 0.1114027E+00, 0.1095039E+00 /

  data h / &
        0.3920617E-01, 0.3932705E-01, 0.3950999E-01, 0.3975703E-01, &
        0.4007093E-01, 0.4045533E-01, 0.4091481E-01, 0.4145507E-01, &
        0.4208311E-01, 0.4280748E-01, 0.4363863E-01, 0.4458932E-01, &
        0.4567523E-01, 0.4691571E-01, 0.4833487E-01, 0.4996298E-01, &
        0.5183859E-01, 0.5401138E-01, 0.5654656E-01, 0.5953130E-01, &
        0.6308489E-01, 0.6737503E-01, 0.7264544E-01, 0.7926471E-01, &
        0.8781922E-01, 0.9930398E-01, 0.1155599E+00, 0.1404344E+00, &
        0.1836142E+00, 0.2790016E+00, 0.7010474E+00 /

  data t / &
        0.7673828E-03, 0.2306870E-02, 0.3860618E-02, 0.5438454E-02, &
        0.7050699E-02, 0.8708396E-02, 0.1042357E-01, 0.1220953E-01, &
        0.1408125E-01, 0.1605579E-01, 0.1815290E-01, 0.2039573E-01, &
        0.2281177E-01, 0.2543407E-01, 0.2830296E-01, 0.3146822E-01, &
        0.3499233E-01, 0.3895483E-01, 0.4345878E-01, 0.4864035E-01, &
        0.5468334E-01, 0.6184222E-01, 0.7047983E-01, 0.8113195E-01, &
        0.9462444E-01, 0.1123001E+00, 0.1364980E+00, 0.1716886E+00, &
        0.2276241E+00, 0.3304980E+00, 0.5847031E+00 /

  CALL RANDOM_NUMBER(u)
  if ( u <= 0.5E+00 ) then
    s = 0.0E+00
  else
    s = 1.0E+00
  end if
  u = 2.0E+00 * u - s
  u = 32.0E+00 * u
  i = int ( u )
  if ( i == 32 ) then
    i = 31
  end if
!
!  Center
!
  if ( i /= 0 ) then

    ustar = u - real ( i )
    aa = a(i)

    do

      if ( t(i) < ustar ) then

        w = ( ustar - t(i) ) * h(i)

        y = aa + w

        if ( s /= 1.0E+00 ) then
          snorm = y
        else
          snorm = -y
        end if

        return

      end if

      CALL RANDOM_NUMBER(u)
      w = u * ( a(i+1) - aa )
      tt = ( 0.5E+00 * w + aa ) * w

      do

        if ( tt < ustar ) then
          y = aa + w
          if ( s /= 1.0E+00 ) then
            snorm = y
          else
            snorm = -y
          end if
          return
        end if

        CALL RANDOM_NUMBER(u)

        if ( ustar < u ) then
          exit
        end if

        tt = u
        CALL RANDOM_NUMBER(ustar)

      end do

      CALL RANDOM_NUMBER(ustar)

    end do
!
!  Tail
!
  else

    i = 6
    aa = a(32)

    do

      u = u + u

      if ( 1.0E+00 <= u ) then
        exit
      end if

      aa = aa + d(i)
      i = i + 1

    end do

    u = u - 1.0E+00
    w = u * d(i)
    tt = ( 0.5E+00 * w + aa ) * w

    do

      CALL RANDOM_NUMBER(ustar)

      if ( tt < ustar ) then
        y = aa + w
        if ( s /= 1.0E+00 ) then
          snorm = y
        else
          snorm = -y
        end if
        return
      end if

      CALL RANDOM_NUMBER(u)

      if ( u <= ustar ) then
        tt = u
      else
        CALL RANDOM_NUMBER(u)
        w = u * d(i)
        tt = ( 0.5E+00 * w + aa ) * w
      end if

    end do

  end if

end


function sexpo ( )

!*****************************************************************************80
!
!! SEXPO samples the standard exponential distribution.
!
!  Discussion:
!
!   This procedure corresponds to algorithm SA in the reference.
!
!  Licensing:
!
!    This code is distributed under the GNU LGPL license.
!
!  Modified:
!
!    31 March 2013
!
!  Author:
!
!    Original FORTRAN77 version by Barry Brown, James Lovato.
!    FORTRAN90 version by John Burkardt.
!
!  Reference:
!
!    Joachim Ahrens, Ulrich Dieter,
!    Computer Methods for Sampling From the
!    Exponential and Normal Distributions,
!    Communications of the ACM,
!    Volume 15, Number 10, October 1972, pages 873-882.
!
!  Parameters:
!
!    Output, real ( kind = 4 ) SEXPO, a random deviate from the standard
!    exponential distribution.
!
  implicit none

  real ( kind = 4 ) a
  integer ( kind = 4 ) i
  real ( kind = 4 ) q(8)
  real ( kind = 4 ) r
  real ( kind = 4 ) sexpo
  real ( kind = 4 ) u
  real ( kind = 4 ) umin
  real ( kind = 4 ) ustar

  save q

  data q / &
       0.6931472E+00, &
       0.9333737E+00, &
       0.9888778E+00, &
       0.9984959E+00, &
       0.9998293E+00, &
       0.9999833E+00, &
       0.9999986E+00, &
       0.9999999E+00 /

  a = 0.0E+00
  CALL RANDOM_NUMBER(u)

  do

    u = u + u

    if ( 1.0E+00 < u ) then
      exit
    end if

    a = a + q(1)

  end do

  u = u - 1.0E+00

  if ( u <= q(1) ) then
    sexpo = a + u
    return
  end if

  i = 1
  CALL RANDOM_NUMBER(ustar)
  umin = ustar

  do

    CALL RANDOM_NUMBER(ustar)
    umin = min ( umin, ustar )
    i = i + 1

    if ( u <= q(i) ) then
      exit
    end if

  end do

  sexpo = a + umin * q(1)

  return
end




